function [Np,P,phi,PTC,outsigma,R,Vb,Vbsaveout,zeta] = dreumWTP(chunksize,H,K,Np,Nu,quota,...
    rho,insigma,tau,tcidx,tcSplit,theta0,tolInner,a,Vbsavein)
% Objective function for DREUM counterfactual policy analysis. This function uses 
% the shares of applicants with p preference points to hunt j (insigma) as an 
% input then uses this to estimate the predicted shares of applicants (outsigma).

Ncpu = 200;
alpha = theta0(end-1);
s = theta0(end);
pit = [theta0(4:5) (1 - sum(theta0(4:5)))];

% Calculate probability of success for each hunt/round
% [phi,PMarkov] = dreumProbWTP(H,K,Np,quota,insigma);
[phi,PMarkov] = dreumProbWTP2(H,K,Np,quota,insigma);
impossible = phi == 0; % Identify hunt/pp combinations for which \phi_jp = 0


% Calculate innner nest, choice probabilities
Pout = cell(1,Ncpu);
zetaout = cell(1,Ncpu);
Rout = cell(1,Ncpu);
Vbout = cell(1,Ncpu);
Vbsaveout = cell(1,Ncpu);
parfor idx1 = 1:Ncpu

    tctp=tcSplit{idx1};
	temp=Vbsavein{idx1};
	Vbsavetp=cell(1,tau);
	for idx2 = 1:tau 
		Vbsavetp{idx2} = temp(:,:,:,idx2);
    end
        
    % Calculate maximized net present value from bear hunting 
%     Vb = dreumInnerWTP(chi0,eta,H,impossible(2:end,:),K,mu,Nu,phi,PMarkov,rho,...
%         tau,tc,tolInner);
    [Vb,Rout{idx1}] = dreumInner5(H,impossible,K,tau,phi,PMarkov,rho,tctp,theta0,tolInner,Vbsavetp);
    Vbout{idx1} = Vb;
    temp = zeros(H+1,K,size(tctp,1),tau);
	for idx2 = 1:tau
		temp(:,:,:,idx2) = Vb{idx2};
	end
	Vbsaveout{idx1} = temp;
    
    % Calculate type probabilities (pi_taup) and site chioce probabilities (P) 
%     pi = [1 - pi(1) - pi(2), pi(1), pi(2)];
%     [~,~,P,PTC,Pp,zeta] = dreumStationary_6(alpha,data,H,impossible,K,Nu,phi,...
%         pi,tau,tcidx,Vb);
    [Pout{idx1},zetaout{idx1}] = dreumStationary_6_WTP(alpha,H,impossible,K,phi,s,tau,Vb);
    
%     fprintf('CPU:\t%d\n',idx1)

end

% save('Vbsave.mat','Vbsave')

% Combine parallel outputs into a single matrix for each variable
P = zeros(H+1,K,Nu,tau);
zeta = zeros(K,1,Nu,tau);
Vb = zeros(H+1,K,Nu,tau);
R = zeros(K,1,Nu,tau);
Nuctr2 = 0;
for idx1 = 1:Ncpu
    n = chunksize(idx1);
    for idx3 = 1:n 
        Nuctr2 = Nuctr2+1;
        for idx2 = 1:tau
            P(:,:,Nuctr2,idx2) = Pout{idx1}{idx2}(:,:,idx3);
            zeta(:,:,Nuctr2,idx2) = zetaout{idx1}{idx2}(:,:,idx3);
            Vb(:,:,Nuctr2,idx2) = Vbout{idx1}{idx2}(:,:,idx3);
            R(:,:,Nuctr2,idx2) = Rout{idx1}(:,idx3,idx2);
        end

    end

end


% Calculate new # of applicants w/ p preference points
PTC = zeros(Nu,1);
for idx1 = 1:Nu
    PTC(idx1) = sum(tcidx==idx1)/size(tcidx,1);
end

PpTC = zeros(Nu,K);
for idx2 = 1:Nu
    for idx3 = 1:K
        PpTC(idx2,idx3) = zeta(idx3,1,idx2,1)*pit(1) ...
            + zeta(idx3,1,idx2,2)*pit(2) + zeta(idx3,1,idx2,3)*pit(3);
    end
end

Pp = PpTC'*PTC;
Np = a.*Np(1,1:K)+(1-a).*(Pp').*55454; 


% Randomize over impossible sites
impCt = impossible./repmat(sum(impossible),[H+1,1]);
for idx1 = 1:tau
    for idx2 = 1:Nu
        y = P(:,:,idx2,idx1);
        y(impossible) = y(impossible).*impCt(impossible);
        P(:,:,idx2,idx1) = y;
    end 
end


% Calculate estimated shares of applicants with k preference points that visit each site
x = zeros(H,K,Nu,tau);
for idx1 = 1:tau
    for idx2 = 1:H
        for idx3 = 1:K
            for idx4 = 1:Nu
                x(idx2,idx3,idx4,idx1) = P(idx2+1,idx3,idx4,idx1)...
                    *zeta(idx3,1,idx4,idx1)*pit(idx1)*PTC(idx4)/Pp(idx3);
            end
        end
    end
end

x = sum(x,4);
% x = x{1} + x{2} + x{3};
outsigma = sum(x,3);
outsigma = [1 - sum(outsigma); outsigma];